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We investigate cusp singularities in f{R) gravity, especially for Starobinsky and Hu-Sawicki dark 
energy models. We illustrate that, by using double-null numerical simulations, a cusp singularity can 
be triggered by gravitational collapses. This singularity can be cured by adding a quadratic term, 
but this causes a Ricci scalar bump that can be observed by an observer outside the event horizon. 
Comparing with cosmological parameters, it seems that it would be difficult to see super-Planckian 
effects by astrophysical experiments. On the other hand, at once there exists a cusp singularity, it 
can be a mechanism to realize a horizon scale curvature singularity that can be interpreted by a 
firewall. 
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I. INTRODUCTION 


It is generally believed that Einstein gravity is not a complete theory and therefore it should be modified. In 
this regard, the introduction of higher curvature terms has been a natural choice by many researchers, motivated 
partly by its connection with fundamental concepts such as string theory |l[, 2 . Among various contributions 
of higher curvature terms, a function of the Ricci scalar is the most simplest extension; such model can be 
transformed to a scalar-tensor gravity model as well as Einstein gravity via a conformal transformation fefl. 
This is commonly referred to as f(R) gravity. 

Since f{R) gravity is transformable to a Brans-Dicke type theory Q, f{R) gravity is in general easier to 
control compared to the other general higher curvature corrections. This means that the behavior of the Ricci 
scalar can be understood from the dynamics of the effective gravitational coupling, i.e., the Brans-Dicke field. 
Thus, f(R) gravity can provide a natural motivation for the choice of certain shape of a scalar field potential 
in the Einstein frame that may be useful for an inflationary model Q or a dark energy model for late time 
cosmology j|| 9]. 

The dynamics of the Brans-Dicke field will be determined by its potential that depends on the choice of /(/?). 
There is, however, a shortcoming of such an approach. If the shape of f(R) of a model is unwisely chosen, then 
it may result in a cusp in the potential of the Brans-Dicke field; this may imply the divergence of the Ricci 
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scalar [l^]. If a region has a divergent Ricci scalar, then it should mean a singularity. However, since this is 
due to an unwise choice of f(R ), such singularity may be fictitious and it may mean that there should be a 
way to regularize it. We call this as the f(R )~induced singularity or the cusp singularity [id). 

The follows are known results in the literature: 

- Many dark energy models of f(R) gravity contain a cusp singularity. In order to obtain a stable and 
non-vanishing (de Sitter) local minimum, one often has to fabricate the shape of the effective potential 


of the Brans-Dicke field; equivalently, one has to tune a form of f(R) |8j, |9|. For this, the R n (n > 0) 
corrections have often been introduced, which nevertheless become the origin of the trouble. 

If the correction term is quadratic, i.e., ~ R 2 , then this cusp singularity can be ameliorated fic| . Such a 
correction term is highly motivated as a viable model for inflation [g, 7J. 

If the cusp singularity is inevitable, then it may be realized during a gravitational collapse Q Q. 
Furthermore, such a singularity can be located outside the event horizon. Therefore, if it does exist, then 
it can cause a detectable naked singularity. 

From these, a question naturally arises as follows. Even though we cure the cusp singularity, there are still 
possibilities that a gravitational collapse induces a very higher curvature region. Even though it does not 
diverge, if it reaches around the Planck scale, then it will cause Planck-scale physical effects. Is it really 
possible and what will be the parameters that determine such a phenomenology? Regarding this, there can be 
two opinions. 

Cons'. As long as there is a Planck-scale effect, it can result a distinctive outcome to outside the black hole, 
while we do not have any evidences of Planck-scale effects during gravitational collapses. Therefore, we 
should avoid a model that causes a cusp singularity. 

Pros : If the existence of a naked singularity is inevitable, then it implies that there appears Planck-scale effects 
around the horizon scale. This will be useful to understand the unitarity issue of black hole dynamics. 

The discussion on these two opposite opinions is the motivation of this paper. 

In this paper, we study a gravitational collapse of f(R) dark energy models that contain cusp singularities, 
by using numerical simulations. Especially we use double-null formalism [ 13 I to study the formation of charged 
black holes [3, 151 where this can be applied for cosmology [h]] as well as scalar-tensor models [ill 18|. We 
first show that we can cure the cusp singularity problem by adding a ~ R 2 term; in addition, we show that 
there can be a higher curvature region. 

This paper is organized as follows. In SEC. El we illustrate our model as well as its realization to double¬ 
null formalism. In SEC. ED we discuss the details of numerical results and observe phenomena of cured cusp 
singularities. Finally, in SEC. HVl we summarize and discuss observational implications and future applications. 
In this paper, we use the natural units: c = G = h = 1. 
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II. DOUBLE-NULL FORMALISM FOR f(R) GRAVITY 


In this section, we discuss the way to realize double-null formalism for f(R) gravity models. We follow 


previous notations of Q. To realize f(R) gravity to double-null formalism, we especially used a technique 
of 11]. Mathematical details are in Appendix A and its convergence and consistency checks are discussed in 
Appendix B. 


A. ,f{R) gravity 


The action of the f(R) gravity model is presented by 


S= d^X\J—g 


-. n f (7?) T ^matter 

i07T 


(1) 


where R is the Ricci scalar and f(R) is a given function of R. By introducing an auxiliary field, this model can 
be transformed to the Brans-Dicke type model 


s = I d A xv~g m - m )+c 


matter 


( 2 ) 


where V($) = —f(R) + Rf'(R) and $ = /'(i?). This is the u> = 0 limit of the Brans-Dicke theory with a 
potential R(4>). 

In order to induce a gravitational collapse, one convenient way is to introduce a gauge field, A^, that couples 
to a matter field, for which we consider a complex scalar field </>. To implement this, we use the following model 

Q,Q: 


-S'bd = / d A Xy/^g 


1 

167T 




•EM 


where 


£ EM = (</>„ + ieA^) - ieA v J) - -1-F^F^, 


( 3 ) 


( 4 ) 


e is a gauge coupling and F^ v = A„ ;a , — A^ v . We choose w = 0 throughout this paper. 
The Einstein equation is as follows: 

'J'C 

= 8t rT™ + 8 tt^ = 8ttT„„, 

where the Brans-Dicke part of the energy-momentum tensor is 

r l]j.v = (— 9p.v^\pa9 P + ~ 1^9 ip® -a9 P ^ — 9pv 

and the matter part of the energy-momentum tensor is 


U(£) 

167T 


( 5 ) 


( 6 ) 


Tfiu — n A </> ;M <?V) T n ( <t>;iii^A v ^> + c^.jjieA^cj) -f- < p-^ieA^cj)^ 


+ — F^pF v p + e A^A^cjxj) + C g^. 
47T 


( 7 ) 
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Field equations for <I>, <j>, and A M are as follows: 

8ir 


1 


<T> flM"_7"C _ 

;fM/y 3 + 2 u 3 + 2 

4>-,LU'9 fJ ' 1 ' + ( 2 0 ;m + ie-A^cj)) + ieA^g^cf) = 0 


2V 1 = 0, 


■^-F l '^v-ie(j)(<p. li -ieA l ^)+ie^{(j). tl + ieA t j lj (jj) = 0 , 

where T c = T c ^ . For convenience, we define the effective potential [/(4>) by 

/ <E> 

(W'(3-) - 21/(3-)) d3>, 

where this determines the dynamics of the field. 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 


B. Initial conditions and free parameters 


We use the double-null coordinate 


ds 2 = — a 2 (u, v)dudv + r 2 (u, v)dfl 2 


( 12 ) 


for numerical simulations. Here, u is the retarded time, v is the advanced time, and dll 2 is the two-sphere. For 
numerical integration, it is convenient to change all equations by a set of first order differential equations. For 
this purpose, we introduce the following notations [llj,llij, lla, llal : the metric function a, the radial function r, 
the Brans-Dicke field $, and a scalar field s = \/47r</>, and define 

OL OL 

h=—, d = —, f = r u , g = r tV , W=$ >u , Z = w = s, u , z = s >v . (13) 

a a 

In addition, in order to avoid the complicated relation between f(R) and l/(<f"), we use the evolution of Ricci 
scalars (see Appendix A). Then we need initial conditions for all functions (a, h, d, r, /, g , $, W, Z, s, w, z, a, q , R ) 
at the initial u = u\ and v — Vi surfaces, where we set iq = Vi = 0. 

We have a gauge freedom to choose the initial r function, i.e., we have a freedom to choose distances 
between null lines. We choose r(0, 0) = ro, f(u,0) = r u o, and g(0,v) = r v o, where r u o < 0 and r v o > 0 
in order to implement that r increases for an out-going direction and decreases for an in-going direction. 
We assume that $ is asymptotically a constant 4>o at a local minimum of the effective potential, so that 
$cW'($o) - 2+($o) = 2f(Ro) - Rof(Ro) = «• Then, $(u, 0) = $(0,u) = $ 0 , W(u, 0) = Z(0 ,w) = 0, and 
R(u,0) = R(0,v) = Rq. In addition, we need the following information. 

In-going null direction: We choose s(u, 0) = 0 and w(u,0) = h(u, 0) = a(u,0) = q(a, 0) = 0. Since the 
asymptotic Misner-Sharp mass function 

1/ 

r“ 6$ 


m(u,v) = - ( l+4^“ r 
2 \ ol 


w ,v q y 2 
2 h 12 — ~ r 


(14) 


should vanish at u = v = 0, it is convenient to choose r u o = —1/2 and r v o = 1/2, and 


a (0,0) = 1- 


V($o) 

6$ 0 


- 1/2 


(15) 


We need more information to determine d,g, Z , and z at the v = 0 surface. We obtain d from Eq. S3, 
g from Eq. (1481) . Z from Eq. (1491) . and z from Eq. (ED. 
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FIG. 1: V($) and [/($) for f s (R). 




FIG. 2: V($) and !/($) for /hs(A). 


Out-going null dierction: In order to derive a gravitational collapse, we need to choose a shape of an in¬ 
going scalar field s(0, v). In principle, we can choose this scalar field profile arbitrarily, but for convenience, 
we use 


s(zij, v) = A sin 2 I 7r 


V — Vi 
V{ — Vi 


cos 2ir 


V — V[ 


+ i cos ( 2-7T 


V — Vi 


, v . ( 16 ) 

Vf — Vi J \ Vi — Vi 

for 0 < v < V{ and s(0,i>) = 0 otherwise, where vt is the width of the pulse, A is the amplitude, and 5 is a 
free parameter to tune the phase of the complex scalar field 11, 141. Taking the derivative with respect 
to v. we obtain z(ui,v). From Eq. (l46l) . we find d = r\z\ 2 /2g<& on the u = 0 surface. By integrating d 
along v, we also get a(0,v). However, we need more information for h, f,W,w,a, and q at the u = 0 
surface. As shown in Appendix A, we obtain h from Eq. m > / fr° m Eq- (BHD - W from Eq. (|i5]) , w from 
Eq. (l55l) . a from Eq. (l53l) . and q from Eq. (l5il) . 

Finally, we choose r 0 = 10, 6 = 0.57T, uj = 0, Vf = 20, and leaving A and e free. 

C. Choice of f(R ): cusp singularities and their remedies 


In addition to the initial conditions, we need to specify the detailed models of f(R). In this paper, we are 
especially interested in four models: the Starobinsky (dark energy) model jg], the Hu-Sawicki model [9I, and 
their cured forms by adding ~ R 2 corrections. 
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a. Starobinsky model The Starobinsky (dark energy) model is as follows 


fs{R) = R + XR 0 




( 17 ) 


where A, i?o, and n are free parameters. We choose A = 2, i?o = 0.001, and n = 2 (FIG. [[]). We can see that 
there is a stable equilibrium around $o — 0.9932 and Vo ~ 0.001966. 

In the large R limit, by using the Taylor expansion, we observe 


fs(R) - R-XR 0 + XR 2 0 n+1 ^ + :., 

fs(R) * 1 - 2 « Ai? o” +1 ^TT + -, 
fs( R ) - 2n(2n + l)A<+ 1 -J^ + ..., 

lim £/($) = lim j — 2V A (4>)^ 

= lirn^ j R ( 2 f s (R) - Rfs(R)) fg(R)dR 
r°° dR 

“ J i? 2 ”+ 1 ' 


(18) 

(19) 

( 20 ) 
( 21 ) 
( 22 ) 
(23) 


If n > 0, then the effective potential [/(4>) of the $ = 1 limit is finite, while the effective force term 2/ — Rf = 
— 2V (hence, the limit where the Ricci scalar diverges) is infinite. This shows the possibility to reach the 
cusp singularity. 

b. Hu-Sawicki model We second consider the Hu-Sawicki model 


/hs(-R) = R ~ X 


C| R r ' 


(24) 


C 2 R n + X r ‘ 

We choose A = 1, C\ = 0.01, C 2 = 40, and n = 1 (FIG. [2J. We can see that there is a stable equilibrium 
around 4> 0 ~ 0.99 and Vo — 0.001. 

In the large R limit, by dehning A = A /C^ n , we observe 

ACi ACi A" 


/hs(-R) — R — 
fks(R) — 1 — n 


/hs(-R) - n(n+ 1) 


C 2 C 2 R n 
XC\ A” 

C 2 R^ 1 + ' 
ACi A 71 


C 2 R n + 2 


Jim [/(<!>) oc J 


dR 

R n+l ' 


(25) 

(26) 

(27) 

(28) 


Therefore, as was in the Starobinsky model, if n > 0, there appears a cusp singularity. 

c. Cure of the cusp singularity by inserting (1/2 )cR 2 In order to avoid the f(R)- induced singularity and 
to study the general feature of the f(R) gravity, we add a correction term that can cure the cusp singularity 














id ]. Here, we study the following ‘cured’ models 1 : 


/cure.s(-R) = k(R)+\cR 2 , (29) 

/cure.Hs(-R) = /hS {R) + ^ cR 2 , (30) 

where c is an arbitrary constant. If R is large enough, then 

f(R)~R+^cR 2 . (31) 

Since /'(?/>) ~ 1 + ctj) = <b, we obtain 

V($) - J7($) — ^ (^ - l) 2 , (32) 

and hence there is no cusp singularity now. For stability, we need to choose c > 0. Throughout this paper, we 
leave c as a constant. 


III. RESULTS AND INTERPRETATIONS 

In this section, we first report numerical results on cusp singularities and their regularizations. Second, we 
discuss the pros and cons of their existence and the possibility of their observations. 

A. Numerical results 

d. Numerical results FIG. [3] shows the cases of naked singularities in f(R) dark energy models 2 (upper 
left is for fs{R) and lower left is for fus{R))- However, after one includes the R 2 regularization term, one can 
extend beyond the cusp singularity (upper right is for f C ure,s(R) and lower right is for / cur e,Hs(-R))- Before the 
singularity, the causal structures of f culet s{R) and f cule ,Hs(R ) are almost the same as the causal structures of 
fs{R) and /hs(-R). However, around the would-be cusp singularity, the causal structures show some difference 
between original and regularized models. One typical behavior of the regularized causal structure is that, soon 
after the would-be cusp singularity, the apparent horizon begins to oscillate. This is due to the dynamics of the 
Brans-Dicke field: FIG. [J] explicitly shows $ for the cured cases. Without the R 2 regularization term, we would 
have $ = 1 at R = oo. On the other hand, after we add the R 2 term, $ = 1 is no more a curvature singularity. 
One can clearly see that the Brans-Dicke field becomes larger and larger around the would-be cusp singularity. 
Due to the large perturbations of the Brans-Dicke held, as it goes to the local minimum, there are oscillatory 
effects of the Held (as well as energy-momentum tensors) and this makes the apparent horizon oscillatory. 


1 Here, we abuse the parameter c as a constant. 

2 In this paper, we used <5 = 0.57T, which this is the case when the charge-mass ratio of a black hole ( Q/M) is maximum EH- 
There is a tendency that as the charge-mass ratio decreases, the effect of the cusp singularity decreases and the cusp singularity 
may appear inside the event horizon. For details, see [ll| . 
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A=0.25, e=0.3, fs(R) A=0.25, e=0.3, MR) 



FIG. 3: Upper: causal structures for fs{R) and f C uie,s{R) with A = 0.25, e = 0.3, and c = 0.001. Lower: causal 
structures for /hs(A) and / cur e,Hs(A) with A = 0.25, e = 0.6, and c = 0.001. 



FIG. 4: <f> for f C uie,s{R) with (A = 0.25, e = 0.3, c = 0.001) and f C ure,Hs(R) with (A = 0.25, e = 0.6, c = 0.001). 


e. Ricci scalar bump Another salient feature is the behavior of the Ricci scalar. If there is no R 2 correction 
term, then the Ricci scalar reaches infinity. However, due to the R 2 correction term, the Ricci scalar becomes 
bounded. The point of FIG. [5] is that, even though the Ricci scalar is bounded, it can increase to a sufficiently 
large value. The blue colored strips in FIG. [5] show such behaviors. The Ricci scalar can become larger than 
the Planck scale (~ 1 in our setup), which can in principle be detected by an observer outside the black hole. 

We can further investigate the size of the Ricci scalar bump by varying the parameter c. Although there is no 
canonical way to define the size of the bump, we may appreciate the maximum value of the Ricci scalar through 
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A=0.25, e=0.3, fs(R) 



A=0.25, e=0.3, £„„,s(R) A=0.25, c=0.6, fHs(R) A=0.25, e=0.6, W.hs(R) 



FIG. 5: Left: the causal structure for fs{R) with (A = 0.25, e = 0.3) and the Ricci scalar R for /cure,s (A) with 
{A = 0.25, e = 0.3, c = 0.001). Right: the causal structure for /hs(A) with (A = 0.25, e = 0.6) and the Ricci scalar R 
for /cure,Hs(R) with (A = 0.25, e = 0.6, c = 0.001). In these figures, in the white-colored region, the Ricci scalar is close 
to zero, while it becomes relatively larger values around the blue-colored region. 




V 


V 


FIG. 6: Ricci scalar on the event horizon (as functions of v) for f C ure,s(R) (left) and / CU re,Hs(A) (right), by varying 
c = 0.001, 0.005, 0.01, 0.1. 

the event horizon 3 . Even though we vary c, its back-reaction to the asymptotic mass remains unchanged, and 
hence the location of the event horizon is not appreciably shifted (though the apparent horizon is sensitively 
changed by varying c). FIG. [6] summarizes our results: as c increases, the maximum value of the Ricci scalar 
decreases, where the dependence is proportional to ~ 1 /y/c (FIG. [3). 

/. Causal structures FIG. [S] is the summary of the causal structure of gravitational collapses in f(R) 
gravity. (A) shows the case when the cusp singularity is inside the event horizon. (B) is the case when the 
cusp singularity is outside the event horizon. If one inserts the R 2 regularization term, then one can extend 


3 In principle, the event horizon is not possible to determine unless the future infinity is known. However, the apparent horizon 
approaches to a null direction and this will be the event horizon if there is no incoming/outgoing matter around the horizon. In 
this regard, we determine the location of the event horizon from the asymptotic null direction of the apparent horizon. 
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FIG. 7: The linear dependence between l/\fc and i? max on the event horizon. 


(A) 


(B) 


(C) 




FIG. 8: (A) f(R )~induced cusp singularity is inside the event horizon [llj. (B) As the mass-charge ratio increases, 
the cusp singularity can be outside the event horizon. (C) If a ~ R 2 term is added, then the cusp singularity can be 
removed but a Ricci scalar bump would appear. This can be located outside the horizon and hence can be observed by 
the outside observer (blue-colored region). 


the causal structure as denoted in (C). However, due to the Ricci scalar bump, the high curvature (perhaps, 
super-Planckian) effects can in principle be observed outside the event horizon (in the blue shaded region). 

B. Cons: can we see super-Planckian effects for realistic cosmological models? 

In this subsection, let us try a crude estimation of the Ricci scalar bump. If [/'($) = 0, then there appears a 
stable local minimum at <f>o such that R( ( I > o) = A will be the value of an effective cosmological constant. Then, 
there exists a typical background Ricci scalar R 0 ~ A. Around R ~ Rq, the f(R) dark energy model (before 
cured) is a good description. On the other hand, due to the perturbation of the field during the gravitational 
collapses, R can increase and if it reaches a certain bound, then R+ ( c/2)R 2 becomes a better approximation. 
At once the theory is approximated by the R + ( c/2)R 2 model, which effectively restricts the increase of the 
Ricci scalar. Therefore approximately, the maximum possible Ricci scalar R max occurs when the correction 
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FIG. 9: Conceptual diagram of the effective potential U(R). The red curve corresponds to the f(R) dark energy models 
with the cusp singularity; one can reach R = oo with a finite energy. The green curve corresponds to the R + ( c/2)R 2 
model. The blue curve is the cured model via the R 2 term, where it is similar as the red curve for the small curvature 
region and the green curve for the large curvature region. Between two regions, there is an effective boundary which is 
on the order of R ~ (A/c) 1 ^ 2 . 


term (c/2)i? 2 becomes comparable to the leading term R, or A ~ ci?^ ax . In FIG. [3l A ~ R 0 ~ 10 -3 and 
c ~ 10 -3 . Therefore, R max ~ 0(1) and this is shown by FIG. [5] 

FIG. [9] gives a conceptual picture. The R 2 term can cure the cusp singularity and hence we obtain the 
blue colored curve, where it follows the red colored curve around the small curvature regime and follows the 
green colored curve for the large curvature regime. Basically, R ~ (A/c) 1 / 2 becomes the boundary between two 
different regimes. The gravitational collapse perturb the scalar field away from the local minimum R -C (A/c) 1 / 2 
to become R ~ (A/c) 1 / 2 . 

Now the question is, what are the proper values of c? Note that the effective mass of the scalar field is 
to 2 = 1/c. Therefore, there is a maximum bound of c such that c < m,,. where Wbd is the experimental bound 
of the scalar field mass. Note that this condition can be weakened by the Chameleon mechanism 191. As we 

n 

also want to include the inflation scenario in the f(R) gravity [6|, we find that the observational constraints 

Although this is not very conclusive, a reasonable interpretation is 
~ 10 -66 would be way sub-Planckian (see also 


EL 


restrict m ~ ICC 5 in Planck units 
that, since A ~ 10~ 123 -C 1 in terms of the Planck units, R 
discussions in H)- 


C. Pros: horizon as a super-Planckian surface? 

If we regard f(R) gravity as a possible modification of the Einstein gravity, then the cusp singularity appears 
inevitable from Ricci scalar terms. Of course, such behavior is generic and can be found not only for the Ricci 
scalar, but also other types of curvature terms, though these will be technically more difficult to investigate. If 
it is the case, then this means that a gravitational collapse can trigger higher curvature corrections and this can 
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cause a curvature singularity that is not in Einstein gravity. Such a curvature singularity can be located outside 
as well as inside the horizon ((A) and (B) of FIG. EJ). The point is that the areal radius of the singularity is 
non-zero if such a singularity is induced by higher curvature terms. 


Let us focus on the recent discussion on the information loss problem 
complementarity argument 


21 j. According to the black hole 


22j , the unitarity can be preserved for all processes of a black hole, from collapse 


to evaporation, to an asymptotic observer, as well as an in-falling observer based on quantum field theory of 
curved spacetime. However, later it was noticed that black hole complementarity is inconsistent 


hence we need to drop one of the assumptions. The authors in 


23 


24j and 


24j tried to drop general relativity for an 


in-falling observer around the horizon scale, and they said that there exists a firewall around the horizon scale. 
If there exists a firewall and if it can prevent the problem of black hole complementarity, then it should affect 
to outside the event horizon [25 1. 


If there exists a firewall, then it may resolve the tension of black hole complementarity; however, what is the 
mechanism to realize the firewall? Up to now, there is no consensus (for useful comments, see [ 26 , 271). However, 
if we remind the f(R) cusp singularities, then it opens a new interesting possibility. If there exists higher order 
curvature corrections from quantum gravitational effects, such effects can be triggered by gravitational collapses. 
Note that if it does not have a cusp singularity, then it will not be triggered by gravitational collapses, since 
gravitational collapses are in general sub-Planckian. 

If there is a horizon scale curvature singularity that prevents an in-falling observer, can we identify this as 
a firewall? Of course, we do not have enough evidences whether a cusp singularity should appear for every 
cases. However, at least, this is one of logical possibilities; also, the properties (causal structures) of the cusp 
singularity induced by f(R) seem to be quite consistent with the notion of the firewall. Therefore, it needs 
more investigations for future studies. 


IV. DISCUSSION 


In this paper, we investigated gravitational collapses in f(R) gravity. For some dark energy models, e.g., 
Starobinsky model [§1 and Hu-Sawicki model Q, gravitational collapses can trigger to create a cusp singularity. 
Such a cusp singularity can be cured by adding a quadratic term ~ cR 2 , but still there can exist a Ricci scalar 
bump around the horizon scale. 

If one of f(R) dark energy models is correct, then it is very natural to accept the existence of a cusp 
singularity, or at least, its cured effect as a Ricci scalar bump. Then we can ask this question: can we really 
observe such an effect outside the horizon? 

Regarding this, one may discuss against this: the curvature scale of a Ricci scalar bump will be approximately 
Umax ~ y/A/c, where A is the cosmological constant of the current universe. Since our cosmological constant 
is much smaller than the Planck scale, i.e., Acl, hence R m ax will be much smaller than the Planck scale. 

On the other hand, one may also discuss for this: such a cusp singularity or a Ricci scalar bump can role as a 
modification of general relativity even around the horizon scale. Here, a cusp singularity is important; usually 
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a gravitational collapsing process is sub-Planckian and hence its back-reaction to higher curvature corrections 
will be much weaker. However, if there exists a cusp singularity, then super-Planckian effect can be triggered 
by sub-Planckian gravitational collapses. 


Can we identify this as a firewall 


24] that resolves troubles of black hole complementarity [23J? It is still 


unclear. At once there is a cusp singularity, even though it is inside the horizon, it will be useful to prevent 
inconsistency of black hole complementarity. In addition, the existence of a cusp singularity can be argued that 
it is originated from quantum corrections of gravity. However, the existence of a cusp singularity requires an 
exotic shape of terms, e.g., R ~ n , while their theoretical stability is not clear. 

At least, we open this as a logical possibility. Perhaps, a gravitational collapse can be a good window to 
investigate Planck scale physics. For further justifications or criticisms, we remain as a future work. 
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Appendix A. Implementation of double-null formalism 

The Einstein tensors and energy-momentum tensor components are expressed as follows (ui = 0): 


G u 

G u 

G v 


= -(/,„-2 fh). 


2 r 2 


(4 rf, v + a 2 + 4 fg) 


- — (g,v - 2 gd ), 
r 

Ggs = —4 —- (d, u H- 


rj~i BD _ 
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87r$ 

Z 


rj- tBD _ _ 


(W tU -2hW), 

U gW + fZ a 2 V 
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1 66 
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1 (Z tV -2dZ), 


32tt$’ 


87r4> 

ji 


2tt a 2 $~’“ ++ 

— \ww + ieaCws — ws) + e 2 a 2 ss] , 
47T 

( a }V f 


r 2 V 
1671 -$ ’ 


r r'^ _ _- 

VV 47r" 


rp C 

1 ee 


47ra 2 ’ 

1 


47ra 2 


(wz + zw) + iea{zs — zs ) + 


2 (a, v y 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


After solving a coupled system of equations, we can write all the equations as follows: 


f,u 

/] 'TTV 

= 2 fh - — {W >u - 2 hW) - — T£, 

(45) 

9,v 

T“ 4.7TT 

= 2 gd — 2tp^’ v ~ 2dZ) —T^ v , 

(46) 
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where 


21 = 
58 = 
£ = 


-^-11 ( 9 W + fZ). 


o?V 
8 $ ’ 


Airr 


a 


fg 


i 


C- t - f - « (sH ' + /2) + s* v ’ 


-fZ-gW- 


4irr 


^-T c + (4>W - 2V) 

2 lo7r 


(50) 

(51) 

(52) 


as well as matter field equations 


— 

a 2 q 

2 ^’ 


(53) 

q, v = 

ier 2 

-— (sz-ss), 


(54) 

W,v = 

fz gw iearz 

ieags ie 2 

(55) 

r r r 

r 4 r 2 a qs • 


In addition, in order to avoid a difficulty to convert f(R) to E(<1>), we substitute using the following identities: 


R($) = -f{R) + Rf\R), 

(56) 

^-2V = 2/(i?) - Rf\R). 

(57) 


Therefore, we need an evolution of the Ricci scalar R. To calculate the evolution of R , we further use the 
following identities: 


R.u = 


$ 


R.v = 


4 >. 


(58) 


f'W ’ f"{RY 

Now, equations for a jUV , and s :UV can be represented by a set of first order differential equations. 

We use the same integration method that was used in previous papers 


Runge-Kutta method 


ci at ; lui nia tui 


. We use the second order 


28]. 


Appendix B. Consistency and convergence tests 

In this appendix, we report on the convergence and consistency tests for our simulations. We tested for 
fcuie,s(R) with (-4 = 0.25, e = 0.3, c = 0.001) and f C me,ns{R) with (A = 0.25, e = 0.6 c = 0.001), since these 
cases include all the interesting features of the present paper. 

For consistency, we can check various relations, but the most important non-trivial test for f(R) gravity will 
be the Ricci scalar. We calculated it by using Eq. (j58]l . while the definition of the Ricci scalar is 



We call the latter Ri and the former R 2 , and checked |i?i — R 2 \/\Ri\ around u = 14,15,16. FIG. flOl shows 
that the differences are less than 10 -3 % for the Starobinsky model and 10 _1 % for the Hu-Sawicki model. 

For convergence, we compared finer simulations: 1 x 1, 2 x 2, and 4x4 times finer. In FIG.QT] we see that 
the difference between the lxl and 2x2 times finer cases is 4 times the difference between the 2x2 and 4x4 
times finer cases, and thus our simulation converges to second order. The numerical error is less than ICE 4 % 
for the Starobinsky model and 10 -6 % for the Hu-Sawicki model. 
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FIG. 10: Consistency test for |i?i — _R 2 |/|-Ri| around u = 14,15,16 for f CU re,s(R ) (Left) and / cur e,Hs(R) (Right). 




FIG. 11: Convergence test: |n x i — r 2 X 2 |/r 2X 2 and 4|r2x2 — r- 4 x 4 |/r 4 X 4 around u = 5,10,15 for f C uie,s(R) (Left) and 
/cure,Hs(R) (Right). This shows the second order convergence. 
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